#clean up
rm(list=ls())

#Load Libraries
library(MASS)
library(TSA)
library(lattice)

#Set Working Directory
setwd('/Volumes/MONOGAN/mediaCount/data/')

#source('http://www.utdallas.edu/~pbrandt/code/pests.r')
#Alternate version of Brandt's code:
wald<-function(coef,restrict,covar,t)
 { # Trap non-conformable restrictions */
   if(ncol(restrict)!=nrow(coef))
     {
     stop("WALD TEST ERROR: Restriction and coefficient matrices are non-conformable")
   }
   # Compute Wald statistic and p-value
   waldval<-t(coef)%*%t(restrict)%*%solve(restrict%*%covar%*%t(restrict))%*%restrict%*%coef
   # Small sample correction
   waldval <- ((t-nrow(coef)+nrow(restrict))/t)*waldval
   wp<-1-pchisq(waldval,nrow(restrict));

   list(Wald.stat=waldval,Wald.p=wp);
 }
####################################################################
# parpdgp -- PAR(p) data generation process
# n = number of observations in generated series
# r = AR(p) coefficients
# x = matrix of regressors, including an intercept.
# d = delta, vector of regression coefficients for x.
#

parpdgp<-function(n,r,x=matrix(0,nrow=n),d=as.matrix(0))
  { # Initialize the parameters
    y<-matrix(0,nrow=n)
    r<-as.matrix(r)
    p<-nrow(r)
    xd<-exp(as.matrix(as.matrix(x)%*%as.vector(d)))
    # Loop over DGP
    for(t in 1:n)
#CHANGE
#      { if(t<(p+1))
      { if(t<(p+13))
          { y[t]<-rpois(1,xd[t]) }
        else
#CHANGE
#          { ylag<-t(y[(t-p):(t-1)])%*%rev(r)
          { ylag<-y[t-1]*r[1]+y[t-2]*r[2]+y[t-3]*r[3]+y[t-6]*r[4]+y[t-8]*r[5]+y[t-11]*r[6]+y[t-12]*r[7]
            y[t]<-rpois(1,(ylag + (1-sum(r))*xd[t]))
          }
      }
    return(ts(y,start=1))
  }

####################################################################
# parpfilter -- Extended Kalman filter for PAR(p) model.  This is used
# to compute the filter and the parameters for the log-likelihood
# function, parpllf().

parpfilter<-function(y,x,r,d)
  { # Initialize some parameters
    y<-as.matrix(y)
    r<-as.matrix(r)
    p<-nrow(r)
    n <- length(y)
    xd<-exp(as.matrix(as.matrix(x)%*%as.vector(d)))
    m<-matrix(0,nrow=n)
    s<-m
    # Filter loop
    for(t in 1:n)
#CHANGE
#      { if(t<(p+3))
      { if(t<(15))
#CHANGE
#          { m[t]<-mean(y[1:(p+3)]); s[t]<-var(y[1:(p+3)]);
          { m[t]<-mean(y[1:(15)]); s[t]<-var(y[1:(15)]);
          }
        else
#CHANGE
#          { ylag<-t(y[(t-p):(t-1)])%*%rev(r)
          { ylag<-y[t-1]*r[1]+y[t-2]*r[2]+y[t-3]*r[3]+y[t-6]*r[4]+y[t-8]*r[5]+y[t-11]*r[6]+y[t-12]*r[7]
            m[t]<-ylag + (1-sum(r))*xd[t]
            s[t]<-(sum(r^2)*var(y[1:t-1]) + (1-sum(r^2))*var(xd[1:t-1]))
          }
      }
    return(ts(matrix(cbind(m,s),ncol=2),start=1,names=c("m","s2")))

  }

####################################################################
# parpllf -- PAR(p) log-likelihood function.
# Returns the value of the LLF for the given data and parameters.

parpllf<-function(p,y=y,x=matrix(0,nrow=length(y)))
  { # initialize the parameters
    k<-ncol(x)
    no.param<-length(p)
    r<-as.matrix(p[1:(no.param-k)])
    d<-as.matrix(p[(no.param-k+1):length(p)])
    y<-as.matrix(y); x<-as.matrix(x);
    # Call the filter
    pf<-parpfilter(y,x,r,d);
    # Extract vectors of parameters from filter output
    pf<-pf[(no.param-k+1):nrow(pf),];
    m<-pf[,"m"]; s2<-pf[,"s2"];
    x<-x[(no.param-k+1):nrow(x),];    y<-y[(no.param-k+1):nrow(y)];
    # Compute log-likelihood function
    lp<-(lgamma((s2*m) + y) - lgamma(s2*m) - lgamma(y+1)
         + (s2*m*log(s2)) - ((s2*m + y)*log(1+s2)))
    lp<-sum(lp)
    lp
  }

####################################################################
# Parp -- this is the main function for calling the PAR(p) model
# estimation.

# Inputs
# y = Event count time series you wish to fit.  (can be a matrix or a
# ts object)
# x = matrix of regressors,including an intercept
# formula = R formula format,, typically y ~ x1 + x2 + ...

# parp.init = REQUIRED starting the AR(p) values for the optimizer.  This is
#             the p AR parameters
# init.param = NULL.  This uses the Poisson regression values for
#             starting values unless the user supplies values in the form
#             c(d_1,d_2,...,d_k), where d_1 is the constant.
#

Parp<-function(formula, p=1, parp.init=rep(0.1,p), init.param=NULL, ...)
  {

    model <- model.frame(formula)
    y <- model.extract(model, "response")
    x <- model.matrix(formula)
    k <- ncol(x)
    xnames <- dimnames(x)[[2]]
    colnames(x) <- xnames

    # Do some work to get a vector of starting values.
    if(is.null(init.param))
       {
         init.param <- glm(formula, family=poisson())$coefficients
       }

    # Do a simplex run to get some good starting values
    start.param <- optim(c(parp.init,init.param), parpllf,
                         method="Nelder-Mead", hessian=FALSE,
                         control=list(maxit=800, reltol=0.0001,
                         trace=1, fnscale=-(length(y))),
                         y=y,x=x,...)

    # Estimate the PAR(p) model
    fitted.param<-optim(start.param$par, parpllf, method = "BFGS", hessian=TRUE,
                        control=list(maxit=1000,
                          trace=2, fnscale=-(length(y))),
                        y=y,x=x,...)

    # Retreive the parameters and llf values
    param<-fitted.param$par
    llf <- fitted.param$value

    # Estimate the standard errors and compute z scores
    covar <- -solve(fitted.param$hessian)
    se<-sqrt(diag(covar))
    z<-param/se

    k <- ncol(x)
    k.plus.p <- length(se)
    num.p <- k.plus.p - k
    fitted.par.filter <- parpfilter(y,x,param[1:num.p],param[(num.p+1):k.plus.p])

    residuals <- y - fitted.par.filter[,"m"]
    # Print the results
    coefs <- t(matrix(rbind(param,se,z),nrow=3));
    rownames(coefs) <- c(rep("rho",num.p),colnames(x))
    colnames(coefs) <- c("Parameters","Std. Errors","Z-score")
##     cat(" ","\n")
##     cat("PAR(p) regression output","\n")
##     cat("--------------------------------------------","\n")
##     printCoefmat(coefs)
##     cat("--------------------------------------------","\n")
##     cat("Log-likelihood value  : ", llf, "\n")
    aic <- (-2*llf)+(2*(length(param)-1))
##    cat("AIC                   : ", aic, "\n")
    dof <- length(y) - k.plus.p
##    cat("Degrees of Freedom    : ", dof, "\n")
##    cat("--------------------------------------------","\n")
    R <- diag(1, nrow=p, ncol=k.plus.p)
    wald.test <- wald(matrix(param, ncol=1), R, covar, length(y))
    wald.p <- wald.test$Wald.p
    wald.test <- wald.test$Wald.stat
##     cat("Test for reduction to Poisson regression\n")
##     cat("--------------------------------------------","\n")
##     cat("Wald test for rho_i=0 : ", wald.test, "\n")
##     cat("Wald p-value          : ", wald.p, "\n")
##     cat("--------------------------------------------","\n\n")
    fit <- list(coefs=coefs,
                optim.param=fitted.param,
                param=param,
                covar=covar,
                std.err=se,
                z=z,
                residuals=residuals,
                llf=llf,
                aic=aic,
                dof=dof,
                wald.test=wald.test,
                wald.p=wald.p,
                k=k,
                p=num.p,
                y=y,
                x=x,
                call=match.call(),
                obj.type=c("PAR(p) regression output"))
    class(fit) <- c("pests.output")
    return(fit)
  }


# Print function for the PEWMA and PAR(p) output.

"print.pests.output" <- function(x)
  { cat(" ","\n")
    cat(x$obj.type, "\n")
    cat("\nCall: ", deparse(x$call), "\n\n")
    cat("--------------------------------------------","\n")
    printCoefmat(x$coefs)
    cat("--------------------------------------------","\n")
    cat("Log-likelihood value  : ", x$llf, "\n")
    aic <- (-2*x$llf)+(2*(length(x$param)-1))
    cat("AIC                   : ", x$aic, "\n")
    cat("Degrees of Freedom    : ", x$dof, "\n")
    cat("--------------------------------------------","\n")
    cat("Test for reduction to Poisson\n")
    cat("--------------------------------------------","\n")
    cat("Wald test statistic   : ", x$wald.test, "\n")
    cat("Wald p-value          : ", x$wald.p, "\n")
    cat("--------------------------------------------","\n")

  }

# Holder function for the summary method.
"summary.pests.output" <- print.pests.output

####################################################################
# ACF function standardized for count data...
# Feed in your counts and you can generate some ACFs!

# Computes an ACF for counts by standardizing them.  This is a valid
# method for searching for autocorrelation in counts.  See Cameron and
# Trivedi (Regression Analysis of Count Data, CUP, 1998: 228).

counts.acf <- function(y, lag.max=(length(y)/4),...)
  { z <- y-mean(y)
    ct.acf <- acf(z, lag.max=lag.max,plot=F,...)
    return(ct.acf)
  }


# Box-Ljung Q-text for count data.  Tests for serial correlation ag
# "lag" value.
counts.Q <- function(y, lag=1, ...)
  { z <- y-mean(y)
    test <- Box.test(z,lag=lag,type=c("Ljung-Box"))
    return(test)
  }

#multipliers
parp.multipliers <- function(parp.obj, xvec=matrix(rep(1,parp.obj$k),1,parp.obj$k))
  { p <- parp.obj$p
    k <- parp.obj$k
    ar.coef <- parp.obj$param[1:p]
    x.coef <- matrix(parp.obj$param[(p+1):(p+k)], ncol=1)
    sum.ar.coef <- sum(ar.coef)
    impact.multiplier <- ((1-sum.ar.coef)*exp(xvec%*%x.coef)*x.coef[,1])
    longrun.multiplier <- impact.multiplier/(1-sum.ar.coef)

    # Format
    impacts <- rbind(impact.multiplier, longrun.multiplier)
    colnames(impacts) <- colnames(parp.obj$x)
    rownames(impacts) <- c("Short Run", "Long Run")
    output <- list(impacts=impacts, call=parp.obj$call, changes=xvec)
    class(output) <- c("pestsmultipliers")
    return(output)
  }

# print.pestsmultipliers: printing function for parp.multipliers

"print.pestsmultipliers" <- function(x)
  { print(x$impacts)
  }


########### FLEMMING, BOHTE & WOOD 1997, TABLE 4 ########################
wood.4<-read.csv('woodCHURCH1.csv')

pdf("fbwHist.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("fbwHist.eps",family='ComputerModern',width=6,height=6,pointsize=10)
#histogram(~wood.4$volume,xlab="Church & State (Reader's Guide)")
plot(table(wood.4$volume),type='h',axes=F,xlab="Church & State (Reader's Guide)",ylab="Frequency",xlim=c(0,33))
axis(1,at=seq(0,35,by=5))
axis(2)
abline(h=0,col='gray60')
dev.off()

###Negative Binomial Regression###
wood.table.4.nb<-glm.nb(wood.4$volume~wood.4$stepIllinois+wood.4$stepEngel+wood.4$pulseLynch); summary(wood.table.4.nb)

#correlation between true and fitted
cor(wood.4$volume[-c(1:12)],wood.table.4.nb$fitted.values[-c(1:12)])
cor(wood.4$volume[-c(1:12)],wood.table.4.nb$fitted.values[-c(1:12)])^2

#RMSE
res<-wood.4$volume[-c(1:12)]-wood.table.4.nb$fitted.values[-c(1:12)]
sqrt(sum(res^2))

#Jarque-Bera
jarque.bera.test(wood.table.4.nb$residuals)

#LB
Box.test(wood.table.4.nb$residuals,16,'Ljung-Box')

#long-run effects
100*(exp(wood.table.4.nb$coefficients)-1)

###Identification and PAR Model###
#standardized residuals from NB reg
wood.4.std<-(wood.4$volume-wood.table.4.nb$fitted.values)/sqrt(wood.table.4.nb$fitted.values)

#ACF/PACF
acf(wood.4.std,26)
pacf(wood.4.std,26)

#Fit the PAR model
wood.table.4.parp<-Parp(wood.4$volume~wood.4$stepIllinois+wood.4$stepEngel+wood.4$pulseLynch,p=7); summary(wood.table.4.parp)
wood.4.parp.std<-(wood.table.4.parp$residuals)/sqrt(wood.4$volume-wood.table.4.parp$residuals)
acf(wood.4.parp.std,26)
pacf(wood.4.parp.std,26)
Box.test(wood.4.parp.std,26,'Ljung-Box')
Box.test(wood.4.parp.std,16,'Ljung-Box')

#correlation between true and fitted
par.fit<-wood.4$volume-wood.table.4.parp$residuals
cor(wood.4$volume[-c(1:12)],par.fit[-c(1:12)])
cor(wood.4$volume[-c(1:12)],par.fit[-c(1:12)])^2

#RMSE
sqrt(sum(wood.table.4.parp$residuals[-c(1:12)]^2))

#Jarque-Bera
jarque.bera.test(wood.table.4.parp$residuals)

#LB
Box.test(wood.4.parp.std,16,'Ljung-Box')
Box.test(wood.table.4.parp$residuals,16,'Ljung-Box')

#long-run effects
parp.multipliers(wood.table.4.parp,xvec=matrix(c(1,1,0,0),1,wood.table.4.parp$k))
100*((3.6584103/par.fit[14]))
parp.multipliers(wood.table.4.parp,xvec=matrix(c(1,1,1,0),1,wood.table.4.parp$k))
100*((-.41212476/par.fit[185]))

###Transfer Function###
wood.table.4.tf<-arimax(wood.4$volume, order=c(12,0,0),  xtransf=cbind(wood.4$stepIllinois,wood.4$stepEngel,wood.4$pulseLynch), transfer=list(c(0,0),c(0,0),c(1,0)),fixed=c(NA,NA,NA,0,0,NA,0,NA,0,0,NA,NA,NA,NA,NA,NA,NA)); wood.table.4.tf

wood.table.4.orig<-arimax(wood.4$volume, order=c(0,0,1),  xtransf=cbind(wood.4$stepIllinois,wood.4$stepEngel,wood.4$pulseLynch), transfer=list(c(0,0),c(0,0),c(1,0))); wood.table.4.orig

#correlation between true and fitted
tf.fit<-wood.4$volume-wood.table.4.tf$residuals
cor(wood.4$volume[-c(1:12)],tf.fit[-c(1:12)])
cor(wood.4$volume[-c(1:12)],tf.fit[-c(1:12)])^2

#RMSE
sqrt(sum(wood.table.4.tf$residuals[-c(1:12)]^2))

#Jarque-Bera
jarque.bera.test(wood.table.4.tf$residuals)

#LB
Box.test(wood.table.4.tf$residuals,16,'Ljung-Box')

#long-run effects
100*((3.0806/tf.fit[14]))
100*((6.6354/tf.fit[185]))


###OLS with lags. (Koyck-like.)###
t.volume<-ts(wood.4$volume)
t.Illinois<-ts(wood.4$stepIllinois)
t.Engel<-ts(wood.4$stepEngel)
t.Lynch<-ts(wood.4$pulseLynch)
l.out<-lag(t.volume,-1)
l2.out<-lag(t.volume,-2)
l3.out<-lag(t.volume,-3)
l6.out<-lag(t.volume,-6)
l8.out<-lag(t.volume,-8)
l11.out<-lag(t.volume,-11)
l12.out<-lag(t.volume,-12)

wood.4.b<-ts.union(t.volume,t.Illinois,t.Engel,t.Lynch,l.out,l2.out,l3.out,l6.out,l8.out,l11.out,l12.out)
wood.table.4.koyck<-lm(t.volume~l.out+l2.out+l3.out+l6.out+l8.out+l11.out+l12.out+t.Illinois+t.Engel+t.Lynch, data=wood.4.b); summary(wood.table.4.koyck)

#correlation between true and fitted
cor(wood.4$volume[-c(1:12)],wood.table.4.koyck$fitted.values)
cor(wood.4$volume[-c(1:12)],wood.table.4.koyck$fitted.values)^2

#RMSE
sqrt(sum(wood.table.4.koyck$residuals[-c(1:12)]^2))

#Jarque-Bera
jarque.bera.test(wood.table.4.koyck$residuals)

#LB
Box.test(wood.table.4.koyck$residuals,16,'Ljung-Box')

#long-run effects
1.82371/(1-0.27396-0.11053-0.07072-0.11354-0.10345-0.09201-0.09078)
100*((12.57644/wood.table.4.koyck$fitted.values[14]))
0.09601/(1-0.27396-0.11053-0.07072-0.11354-0.10345-0.09201-0.09078)
100*((.6620923/tf.fit[185]))

###Log-Normal Model###
log.volume<-ts(log(wood.4$volume+.1))
l.log.out<-lag(log.volume,-1)
l2.log.out<-lag(log.volume,-2)
l3.log.out<-lag(log.volume,-3)
l6.log.out<-lag(log.volume,-6)
l8.log.out<-lag(log.volume,-8)
l11.log.out<-lag(log.volume,-11)
l12.log.out<-lag(log.volume,-12)

wood.4.c<-ts.union(log.volume,t.Illinois,t.Engel,t.Lynch,l.log.out,l2.log.out,l3.log.out,l6.log.out,l8.log.out,l11.log.out,l12.log.out)
wood.table.4.lognorm<-lm(log.volume~l.log.out+l2.log.out+l3.log.out+l6.log.out+l8.log.out+l11.log.out+l12.log.out+t.Illinois+t.Engel+t.Lynch, data=wood.4.c); summary(wood.table.4.lognorm)

#correlation between true and fitted
cor(wood.4$volume[-c(1:12)],exp(wood.table.4.lognorm$fitted.values))
cor(wood.4$volume[-c(1:12)],exp(wood.table.4.lognorm$fitted.values))^2

#RMSE
res<-wood.4$volume[-c(1:12)]-exp(wood.table.4.lognorm$fitted.values)
sqrt(sum(res^2))

#Jarque-Bera
jarque.bera.test(wood.table.4.lognorm$residuals)

#LB
Box.test(wood.table.4.lognorm$residuals,16,'Ljung-Box')

#long-run effects
100*(exp(1.582814/(1-0.133729-0.169534-0.084122-0.091406-0.089731-0.105546-0.096142))-1)
100*(exp(-.007266/(1-0.133729-0.169534-0.084122-0.091406-0.089731-0.105546-0.096142))-1)

######### DRAW THE EFFECT OF LYNCH V. DONNELLY ########

##Clarify-like simulations
fbw.4.sim.par<-mvrnorm(n=1000,mu=wood.table.4.parp$coefs[,1],Sigma=wood.table.4.parp$covar,empirical=TRUE)

fbw.4.results.nb<-summary(wood.table.4.nb)
fbw.4.sim.nb<-mvrnorm(n=1000,mu=fbw.4.results.nb$coefficients[,1],Sigma=fbw.4.results.nb$cov.unscaled,empirical=TRUE)

fbw.4.results.koyck<-summary(wood.table.4.koyck)
fbw.4.covar.koyck<-fbw.4.results.koyck$cov.unscaled*(fbw.4.results.koyck$sigma^2)
fbw.4.sim.koyck<-mvrnorm(n=1000,mu=fbw.4.results.koyck$coefficients[,1],Sigma=fbw.4.covar.koyck,empirical=TRUE)

fbw.4.results.lognorm<-summary(wood.table.4.lognorm)
fbw.4.covar.lognorm<-fbw.4.results.lognorm$cov.unscaled*(fbw.4.results.lognorm$sigma^2)
fbw.4.sim.lognorm<-mvrnorm(n=1000,mu=fbw.4.results.lognorm$coefficients[,1],Sigma=fbw.4.covar.lognorm,empirical=TRUE)

fbw.4.sim.tf<-mvrnorm(n=1000,mu=wood.table.4.tf$coef[-c(4,5,7,9,10)],Sigma=wood.table.4.tf$var.coef,empirical=TRUE)

# 'volume' in February 1984=10. This is what all pre-intervention values are set to.

###Transfer function###
#Start with predictions from reported coefficients.
coef<-wood.table.4.tf$coef[-c(1:12,16)]
lynch.initial<-wood.table.4.tf$coef[17]
lynch.decay<-wood.table.4.tf$coef[16]

#Lynch onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
y00<-coef%*%x00
y0<-coef%*%x0
y1<-coef%*%x1
y2<-coef%*%x2+lynch.decay*lynch.initial
y3<-coef%*%x3+(lynch.decay^2)*lynch.initial
y4<-coef%*%x4+(lynch.decay^3)*lynch.initial
y5<-coef%*%x5+(lynch.decay^4)*lynch.initial
y6<-coef%*%x6+(lynch.decay^5)*lynch.initial
y7<-coef%*%x7+(lynch.decay^6)*lynch.initial
y8<-coef%*%x8+(lynch.decay^7)*lynch.initial
y9<-coef%*%x9+(lynch.decay^8)*lynch.initial
y10<-coef%*%x10+(lynch.decay^9)*lynch.initial

#Plot Forecast
Time<-c(-1:10)
mle.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=mle.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-fbw.4.sim.tf[i,-c(1:7,11)]
lynch.initial<-fbw.4.sim.tf[i,12]
lynch.decay<-fbw.4.sim.tf[i,11]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
forecast[i,1]<-y00<-coef%*%x00
forecast[i,2]<-y0<-coef%*%x0
forecast[i,3]<-y1<-coef%*%x1
forecast[i,4]<-y2<-coef%*%x2+lynch.decay*lynch.initial
forecast[i,5]<-y3<-coef%*%x3+(lynch.decay^2)*lynch.initial
forecast[i,6]<-y4<-coef%*%x4+(lynch.decay^3)*lynch.initial
forecast[i,7]<-y5<-coef%*%x5+(lynch.decay^4)*lynch.initial
forecast[i,8]<-y6<-coef%*%x6+(lynch.decay^5)*lynch.initial
forecast[i,9]<-y7<-coef%*%x7+(lynch.decay^6)*lynch.initial
forecast[i,10]<-y8<-coef%*%x8+(lynch.decay^7)*lynch.initial
forecast[i,11]<-y9<-coef%*%x9+(lynch.decay^8)*lynch.initial
forecast[i,12]<-y10<-coef%*%x10+(lynch.decay^9)*lynch.initial
}

#Lower and Upper Bounds
mle.p.5<-apply(forecast, 2, quantile,.05)
mle.p.95<-apply(forecast, 2, quantile,.95)

###PAR model###
#Start with predictions from reported coefficients.
rho<-wood.table.4.parp$coefs[1:7,1]
weight<-1-sum(rho)
beta<-wood.table.4.parp$coefs[8:11,1]

#Romer onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
lag00<-c(10,10,10,10,10,10,10)
y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(10,10,10,10,10,10,10)
y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(10,10,10,10,10,10,10)
y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1,10,10,10,10,10,10)
y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2,y1,10,10,10,10,10)
y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3,y2,y1,10,10,10,10)
y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4,y3,y2,10,10,10,10)
y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5,y4,y3,10,10,10,10)
y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6,y5,y4,y1,10,10,10)
y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7,y6,y5,y2,10,10,10)
y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8,y7,y6,y3,y1,10,10)
y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9,y8,y7,y4,y2,10,10)
y10<-weight*exp(beta%*%x10)+rho%*%lag10

#Plot Forecast
Time<-c(-1:10)
par.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=par.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
rho<-fbw.4.sim.par[i,1:7]
weight<-1-sum(rho)
beta<-fbw.4.sim.par[i,8:11]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
lag00<-c(10,10,10,10,10,10,10)
forecast[i,1]<-y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(10,10,10,10,10,10,10)
forecast[i,2]<-y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(10,10,10,10,10,10,10)
forecast[i,3]<-y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1,10,10,10,10,10,10)
forecast[i,4]<-y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2,y1,10,10,10,10,10)
forecast[i,5]<-y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3,y2,y1,10,10,10,10)
forecast[i,6]<-y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4,y3,y2,10,10,10,10)
forecast[i,7]<-y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5,y4,y3,10,10,10,10)
forecast[i,8]<-y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6,y5,y4,y1,10,10,10)
forecast[i,9]<-y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7,y6,y5,y2,10,10,10)
forecast[i,10]<-y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8,y7,y6,y3,y1,10,10)
forecast[i,11]<-y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9,y8,y7,y4,y2,10,10)
forecast[i,12]<-y10<-weight*exp(beta%*%x10)+rho%*%lag10
}

#Lower and Upper Bounds
par.p.5<-apply(forecast, 2, quantile,.05)
par.p.95<-apply(forecast, 2, quantile,.95)


###Negative binomial regression###
#Start with predictions from reported coefficients.
coef<-fbw.4.results.nb$coefficients[,1]

#Romer onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
y00<-exp(coef%*%x00)
y0<-exp(coef%*%x0)
y1<-exp(coef%*%x1)
y2<-exp(coef%*%x2)
y3<-exp(coef%*%x3)
y4<-exp(coef%*%x4)
y5<-exp(coef%*%x5)
y6<-exp(coef%*%x6)
y7<-exp(coef%*%x7)
y8<-exp(coef%*%x8)
y9<-exp(coef%*%x9)
y10<-exp(coef%*%x10)

#Plot Forecast
Time<-c(-1:10)
nb.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=nb.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-fbw.4.sim.nb[i,]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
forecast[i,1]<-y00<-exp(coef%*%x00)
forecast[i,2]<-y0<-exp(coef%*%x0)
forecast[i,3]<-y1<-exp(coef%*%x1)
forecast[i,4]<-y2<-exp(coef%*%x2)
forecast[i,5]<-y3<-exp(coef%*%x3)
forecast[i,6]<-y4<-exp(coef%*%x4)
forecast[i,7]<-y5<-exp(coef%*%x5)
forecast[i,8]<-y6<-exp(coef%*%x6)
forecast[i,9]<-y7<-exp(coef%*%x7)
forecast[i,10]<-y8<-exp(coef%*%x8)
forecast[i,11]<-y9<-exp(coef%*%x9)
forecast[i,12]<-y10<-exp(coef%*%x10)
}

#Lower and Upper Bounds
nb.p.5<-apply(forecast, 2, quantile,.05)
nb.p.95<-apply(forecast, 2, quantile,.95)

###OLS with lags. (Koyck-like.)###
#Start with predictions from reported coefficients.
coef<-fbw.4.results.koyck$coefficients[,1]

#Romer onsets at wave 1. Waves -1-10
x00<-c(1,10,10,10,10,10,10,10,0,0,0)
y00<-coef%*%x00
x0<-c(1,10,10,10,10,10,10,10,0,0,0)
y0<-coef%*%x0
x1<-c(1,10,10,10,10,10,10,10,0,0,1)
y1<-coef%*%x1
x2<-c(1,y1,10,10,10,10,10,10,0,0,0)
y2<-coef%*%x2
x3<-c(1,y2,y1,10,10,10,10,10,0,0,0)
y3<-coef%*%x3
x4<-c(1,y3,y2,y1,10,10,10,10,0,0,0)
y4<-coef%*%x4
x5<-c(1,y4,y3,y2,10,10,10,10,0,0,0)
y5<-coef%*%x5
x6<-c(1,y5,y4,y3,10,10,10,10,0,0,0)
y6<-coef%*%x6
x7<-c(1,y6,y5,y4,y1,10,10,10,0,0,0)
y7<-coef%*%x7
x8<-c(1,y7,y6,y5,y2,10,10,10,0,0,0)
y8<-coef%*%x8
x9<-c(1,y8,y7,y6,y3,y1,10,10,0,0,0)
y9<-coef%*%x9
x10<-c(1,y9,y8,y7,y4,y2,10,10,0,0,0)
y10<-coef%*%x10

#Plot Forecast
Time<-c(-1:10)
ols.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=ols.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-fbw.4.sim.koyck[i,]
x00<-c(1,10,10,10,10,10,10,10,0,0,0)
forecast[i,1]<-y00<-coef%*%x00
x0<-c(1,10,10,10,10,10,10,10,0,0,0)
forecast[i,2]<-y0<-coef%*%x0
x1<-c(1,10,10,10,10,10,10,10,0,0,1)
forecast[i,3]<-y1<-coef%*%x1
x2<-c(1,y1,10,10,10,10,10,10,0,0,0)
forecast[i,4]<-y2<-coef%*%x2
x3<-c(1,y2,y1,10,10,10,10,10,0,0,0)
forecast[i,5]<-y3<-coef%*%x3
x4<-c(1,y3,y2,y1,10,10,10,10,0,0,0)
forecast[i,6]<-y4<-coef%*%x4
x5<-c(1,y4,y3,y2,10,10,10,10,0,0,0)
forecast[i,7]<-y5<-coef%*%x5
x6<-c(1,y5,y4,y3,10,10,10,10,0,0,0)
forecast[i,8]<-y6<-coef%*%x6
x7<-c(1,y6,y5,y4,y1,10,10,10,0,0,0)
forecast[i,9]<-y7<-coef%*%x7
x8<-c(1,y7,y6,y5,y2,10,10,10,0,0,0)
forecast[i,10]<-y8<-coef%*%x8
x9<-c(1,y8,y7,y6,y3,y1,10,10,0,0,0)
forecast[i,11]<-y9<-coef%*%x9
x10<-c(1,y9,y8,y7,y4,y2,10,10,0,0,0)
forecast[i,12]<-y10<-coef%*%x10
}

#Lower and Upper Bounds
ols.p.5<-apply(forecast, 2, quantile,.05)
ols.p.95<-apply(forecast, 2, quantile,.95)

###Log-Normal Model###
#Start with predictions from reported coefficients.
coef<-fbw.4.results.lognorm$coefficients[-c(2:8),1]
rho1<-fbw.4.results.lognorm$coefficients[2,1]
rho2<-fbw.4.results.lognorm$coefficients[3,1]
rho3<-fbw.4.results.lognorm$coefficients[4,1]
rho6<-fbw.4.results.lognorm$coefficients[5,1]
rho8<-fbw.4.results.lognorm$coefficients[6,1]
rho11<-fbw.4.results.lognorm$coefficients[7,1]
rho12<-fbw.4.results.lognorm$coefficients[8,1]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
y00<-exp(coef%*%x00)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y0<-exp(coef%*%x0)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y1<-exp(coef%*%x1)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y2<-exp(coef%*%x2)*(y1^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y3<-exp(coef%*%x3)*(y2^rho1)*(y1^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y4<-exp(coef%*%x4)*(y3^rho1)*(y2^rho2)*(y1^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y5<-exp(coef%*%x5)*(y4^rho1)*(y3^rho2)*(y2^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y6<-exp(coef%*%x6)*(y5^rho1)*(y4^rho2)*(y3^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y7<-exp(coef%*%x7)*(y6^rho1)*(y5^rho2)*(y4^rho3)*(y1^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y8<-exp(coef%*%x8)*(y7^rho1)*(y6^rho2)*(y5^rho3)*(y2^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
y9<-exp(coef%*%x9)*(y8^rho1)*(y7^rho2)*(y6^rho3)*(y3^rho6)*(y1^rho8)*(10^rho11)*(10^rho12)
y10<-exp(coef%*%x10)*(y9^rho1)*(y8^rho2)*(y7^rho3)*(y4^rho6)*(y2^rho8)*(10^rho11)*(10^rho12)

#Plot Forecast
Time<-c(-1:10)
lognorm.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=lognorm.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-fbw.4.sim.lognorm[i,-c(2:8)]
rho1<-fbw.4.sim.lognorm[i,2]
rho2<-fbw.4.sim.lognorm[i,3]
rho3<-fbw.4.sim.lognorm[i,4]
rho6<-fbw.4.sim.lognorm[i,5]
rho8<-fbw.4.sim.lognorm[i,6]
rho11<-fbw.4.sim.lognorm[i,7]
rho12<-fbw.4.sim.lognorm[i,8]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0)
x1<-c(1,0,0,1)
forecast[i,1]<-y00<-exp(coef%*%x00)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,2]<-y0<-exp(coef%*%x0)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,3]<-y1<-exp(coef%*%x1)*(10^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,4]<-y2<-exp(coef%*%x2)*(y1^rho1)*(10^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,5]<-y3<-exp(coef%*%x3)*(y2^rho1)*(y1^rho2)*(10^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,6]<-y4<-exp(coef%*%x4)*(y3^rho1)*(y2^rho2)*(y1^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,7]<-y5<-exp(coef%*%x5)*(y4^rho1)*(y3^rho2)*(y2^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,8]<-y6<-exp(coef%*%x6)*(y5^rho1)*(y4^rho2)*(y3^rho3)*(10^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,9]<-y7<-exp(coef%*%x7)*(y6^rho1)*(y5^rho2)*(y4^rho3)*(y1^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,10]<-y8<-exp(coef%*%x8)*(y7^rho1)*(y6^rho2)*(y5^rho3)*(y2^rho6)*(10^rho8)*(10^rho11)*(10^rho12)
forecast[i,11]<-y9<-exp(coef%*%x9)*(y8^rho1)*(y7^rho2)*(y6^rho3)*(y3^rho6)*(y1^rho8)*(10^rho11)*(10^rho12)
forecast[i,12]<-y10<-exp(coef%*%x10)*(y9^rho1)*(y8^rho2)*(y7^rho3)*(y4^rho6)*(y2^rho8)*(10^rho11)*(10^rho12)
}

#Lower and Upper Bounds
lognorm.p.5<-apply(forecast, 2, quantile,.05)
lognorm.p.95<-apply(forecast, 2, quantile,.95)

###Plot All Forecasts With Error Bands###
pdf("fbwLynch.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("fbwLynch.eps",family='ComputerModern',width=6,height=6,pointsize=10)
Time<-c(-1:10)

#points(y=ols.Forecast,x=Time+.1,pch=23,col='blue')
plot(y=ols.Forecast,x=Time-.26,pch=23,ylim=c(-10,50),xlim=c(-1,11),xlab="Time",ylab="Predicted Count",col='blue', axes=F)
arrows(x0=Time-.26,y0=ols.p.5,x1=Time-.26,y1=ols.p.95, code=3, angle=90, length=.0,lty=2,col='blue')
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 84','Apr 84','June 84','Aug 84','Oct 84','Dec 84'))
box()

#plot(y=nb.Forecast,x=Time-.1,pch=20,ylim=c(0,250),xlim=c(-1,11),xlab="Time",ylab="Predicted Count",col='red', axes=F)
points(y=nb.Forecast,x=Time-.13,pch=20,col='red')
arrows(x0=Time-.13,y0=nb.p.5,x1=Time-.13,y1=nb.p.95, code=3, angle=90, length=.0,lty=3,col='red')

points(y=lognorm.Forecast,x=Time,pch=15,col='purple')
arrows(x0=Time,y0=lognorm.p.5,x1=Time,y1=lognorm.p.95, code=3, angle=90, length=.0,lty=5,col='purple')

points(y=mle.Forecast,x=Time+.13,pch=22,col='forestgreen')
arrows(x0=Time+.13,y0=mle.p.5,x1=Time+.13,y1=mle.p.95, code=3, angle=90, length=.0,lty=4,col='forestgreen')

points(y=par.Forecast,x=Time+.26,pch=24)
arrows(x0=Time+.26,y0=par.p.5,x1=Time+.26,y1=par.p.95, code=3, angle=90, length=.0,lty=1)

legend(x=6,y=40,legend=c('Koyck','Negative Binomial','Log-normal','Transfer Function','PAR'), lty=c(2,3,5,4,1), pch=c(23,20,15,22,24), col=c('blue','red','purple','forestgreen','black'))

abline(h=0, col='gray60')
dev.off()
